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We show, for the first time, how to calculate photonic band structures for metals and other 
dispersive systems using an efficient Order N scheme. The method is applied to two simple periodic 
metallic systems where it gives results in close agreement with calculations made with other tech- 
niques. Further, the approach demonstrates excellent numerical stablity within the limits we give. 
Our new method opens the way for efficient calculations on complex structures containing a whole 
new class of material. 

PACS numbers: 42.70.Qs, 02.70.Bf, 78.20.Bh 

In recent years, there has been much interest in artificially structured dielectrics, otherwise known as photonic 
crystals Q. Structured on the scale of the wavelength of light, these materials promise new and exciting optical 
properties based on the novel dispersion relationships, uj(k), induced by the periodic structure. There is a strong 
analogy here with semiconductor physics, where band gaps in the dispersion relationships for electrons play a key role 
in determining the electronic properties. This analogy has been emphasized by many authors ^|J^]. 
Conceptually, the problem reduces to solving Maxwell's equations, 

VxE=-^ VxH = +^, (1) 
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where the photonic structure may be introduced either through the electric permittivity, D(r, u>) = e(r, cj)E(r, uj), or, 
less commonly, the magnetic permeability, H(r, w) = /i(r, w)B(r, uj). In reciprocal space we have, 



ik x E(k, uj) = +iujB(k, uj), ik x H(k, uj) = -iujT>(k, uj). (2) 
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These equations provide the basis for several computational schemes. Theorists have deployed three main techniques 
for calculation of uj (k): plane wave expansions transfer matrix methods in real space |^||, and, more recently, 

the so called 'Order AT' method of Chan, Yu and Ho fiofl . Currently, the field is dominated by the latter two methods, 
each having its own strengths. As in the electronic case computation times for traditional schemes, such as plane 
wave expansions, scale as N 3 , where N is proportional to the size of the system. The scaling law follows from the 
0(A^ 3 ) operations required to diagonalize an iV x JV matrix. This unfavorable scaling with system size has proved 
very restrictive and there has been an increase in activity to devise new methods which scale more favorably with 
the system size. The optimum scaling possible is clearly O(N) - the time taken to define the problem or to look at 
the answer! Chan etal. showed it is possible to realize this optimal scaling by working in real space and in the time 
domain, provided that the resulting equations are local in space and time. That is to say, the fields at each point in 
space-time can be updated using only the recent values of the fields at points which are close by in real space. In 
the context of photonic band structure calculations, this amounts to requiring that the dielectric function and the 
magnetic permeability are constants independent of k and uj. In fact, Chan's method was not a new idea but had 
been known for some time to the electrical engineering community as the 'finite difference time domain (FDTD)' 
method pd[ | , though it had not previously been applied to the photonic band structure problem. 

However, there is a problem with methods such as FDTD. In the time domain, it is not obvious how to deal with 



a frequency dependent dielectric permittivity, e{uj). For example, nearly- free-electron metals can be modeled by, 

£(«;) = ( 3 ) 

UJ Z 

where uj p is the plasma frequency. This restriction is a major obstacle to progress in the field as Order N methods are 
ideal for treating complex unit cells but are excluded from attacking some of the most interesting systems. Metals, 
for example, show some of the most striking effects when micro-structured. Earlier calculations using the transfer 
matrix method show huge enhancements of local fields in nanostructured arrays of metal cylinders or spheres [ fl5| . 
But despite the clear advantage that an O(N) scaling would give in this calculation, it is impossible to do using a 
standard Order N scheme. 
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The electrical engineering community have known about this problem for some time and progress has been made 
towards solving it [[12 13 1. Here we present a method similar to that of Luebbers [Q to calculate for the first time 
the band structure of various photonic structures containing metallic components in an Order TV scheme. Specifically, 
we treat metallic systems exhibiting a simple plasmon pole, as given in equation ([|). Importantly, this can be done 
without destroying the O(TV) scaling property, and we argue that our methods can be extended to treat more complex 
forms of dispersion. This constitutes a substantial advance to Chan's original Order TV band structure scheme and 
allows us to tackle a whole new class of materials previously excluded from the Order TV method. 

We begin with Maxwell's equations from (0) and approximate them as follows, 



ike (k) x E(k, oj) — +iwn (oj) /z^oH(k, oj), ikh (k) x H(k, oj) — —iwe (w) ££oE(k, oj), 



where, 



k-Ex (k x ) = (+iao) 1 [exp (+ik x a Q ) - 1] , n Hx (k x ) = {-ia ) 1 [exp (-ik x a ) - 1] 

WE ( w ) = (-iry 1 [exp (-iuT) - 1] , w H ( w ) = i+iry 1 [ exp (+iuT) - 1] . 



etc. 



(4) 
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In the case that s and ji are independent of cj, Fourier transformation to (r, t) space yields a set of finite difference 
equations on a simple cubic mesh in real space, of lattice constant oq, and on an interval T in the time domain, as 
discussed in ||f|. The resulting equations allow us to update the H and E fields in the time domain. Starting from, 
say, a random set of fields and the appropriate periodic boundary conditions we iterate for a sufficient number of time 
steps, determined by the frequency resolution desired, and then perform a Fourier transform to obtain the frequencies, 

In this paper, we wish to consider band structures for systems in which the permittivity is a not a constant but a 
function of frequency, 



ike (k) x E(k, oj) — +iiVH (w) /ioH(k, oj), inn (k) x H(k, oj) — —iwe (oj) £ (oj) £oE(k, oj). 
The second equation can be written as 

e -1 (uj) inn (k) x H (k, oj) = —iwE (w) £oE (k, oj) . 
We substitute the single pole metallic form from (^|) and approximate, 
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exp [+i (oj — oj p ) T] 
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exp [i (oj — ojp) nT] 
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exp [+i (oj + ojp) T] 
exp [i (oj + oj p ) nT] 



n=0 



and Fourier transform to give, 



£ E (k, t + T) = £ E (k, t) + iTn H (k) x H (k, t) 
T 2 oj n 



■±[U+{t)-U-{t)] 



where, 

U±(t + T) = n H (k) x H (r, t + T) + exp (±ioj p T) U±(t). 
When transformed into real space, these equations together with 

W)H (k, t + T) = u H (k, t) - iT KE (k) x E (r, t) , 



(8) 



(9) 
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give us the expressions for the frequency dependent Order N method. The third term of equation (|l0|) gives the 
contribution due to the frequency-dependent dielectric constant. It should be noted that this term includes a non- 
local time effect, and is a function of the previous history of the system. 

By expanding equation (^) in powers of oj p T it can be shown that the approximation we have made is good for 



frequencies above about oj 2 T/ vl2. In general, the Order TV method is valid up to frequencies of about coir/a^, with 
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the usual stability critirion that T < ao/^/3co. Together, these define the limits of applicability of our approach. In 
fact, an expansion of equation (||) can be used to produce systematic corrections to the approximation and obtain an 
arbitrary lower frequency limit. In our calculations we found that it was sufficient just to correct for the leading error. 

In general, a dielectric function that depends on uj would give equations which require us to remember all previous 
history of the calculation in order to update the fields, and all computational advantage is lost for the scheme. 
However, for the special case where £^ 1 (w) can be represented by a simple pole, as in equation (||), the fields may be 
rapidly updated by defining U±(t) which themselves may be updated from a knowledge of the immediate past only, 
as shown in ([l0|). This stems from the form of the approximation to the pole which we chose in (||). Clearly, this 
result can easily be generalized to the case where can be represented by a sum over poles. The poles have a 

simple interpretation of internal electromagnetically active modes which can absorb energy. Many inverse dielectric 
functions are well approximated by a sum over poles, particularly where the poles lie off the real axis and therefore 
give rise to rather broad structure in e _1 (w). In so far as dielectric functions can be approximated in this manner, 
we have a general technique for including dispersion in an Order N calculation. 

We now present results for the photonic band structure for both one and two dimensional systems. In Fig. |l|, we 
show the photonic band structure for a one-dimensional system formed of alternating layers of a metal of thickness 
a separated by a vacuum layers of thickness 6, the unit cell for which we divide up into ten equally spaced mesh 
points. The metal was characterized by a single pole plasma as given by equation (^), and we considered propagation 
at normal incidence. We chose lo p cI/2t:co = 1, d = a + 6, , / = 0.1, where w p is the plasma frequency and / is the 
metal filling fraction. A total of 8192 time steps were used in the simulation, with each time step of 0.64 in units of 
m e a^/h. From Fig. |we can see that allowed bands exist in the frequency range lo/oj p < 1 in which the dielectric 
function of the metallic component is negative. The most interesting feature of the band structure is the gap below 
the lowest frequency band which exists for any filling fraction. Also shown in Fig. |^ are results for the same system 
obtained with the transfer matrix method for comparison. The excellent agreement gives confidence that the modified 
Order N method is working correctly. 

In Fig. U we plot the width of the gap below the lowest frequency band as a function of the filling fraction. From 
this figure, we see that the width of this gap increases with the filling fraction. This result agrees with previously 
published results |l|jn|, again confirming the validity of our method. 

For the two dimensional case we calculate the photonic band structures for a system consisting of an infinite array 
of parallel, metal rods of square cross-section, embedded in vacuum. The rods are infinitely long in the z-direction 
and arranged on a square lattice in the x — y plane. The thickness of the rods is a, and their separation <i, and our 
unit cell now consists of 10 x 10 mesh points. For the E and H -polarizations we have E (r, i) — [ 0, 0, £3 (r, t) ] 
and H(r,i) = [ 0, 0, Hs(r,t) ] respectively. In Fig. [|, we present the photonic band structure when the filling 
fraction / = 0.01. Figure [| show the result for the photonic band structure when the filling fraction of the rods is 
/ = 0.04. In both cases we again show transfer matrix results for comparison, which again confirm the accuracy of 
the new method. The key features of these results are the low frequency cut-off in the E-polarized bands and the 
very flat bands at around the plasma frequency in the H-polarization. The cut-off is caused by the collective motion 
of electrons screening the electric field parallel to the rods, below some effective plasma frequency determined by the 
filling fraction. The flat bands are caused by the resonant modes of the individual rods, excited by the electric field 
perpendicular to the rods. The reproduction of these previously observed effects again confirms that our approach is 
capturing the correct physics. 

We have calculated, for the first time using Order N techniques, the photonic band structure for photonic materials 
containing metal. Our work represents a significant extention for Chan's Order N band structure method. We expect 
to be able to handle a wide range of other dispersive materials. Within the stated limits, the scheme exhibits excellent 
numerical stability, in common with Chan's original scheme. Our method retains the advantageous order N scaling 
and is therefore capable of treating complex structures with large unit cells. 
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FIG. 1. Photonic band structure for a periodic series of infinite metallic sheets. The filling fraction / = 0.1, and w p d/2-RC = 1. 
Calculated using the Order N method. 

FIG. 2. Photonic band structure for a periodic series of infinite metallic sheets. The filling fraction / = 0.1, and w p d/2iTC = 1. 
Calculated using the transfer matrix method. 

FIG. 3. The width of the gap below the lowest frequency band as a function of the filling fraction filling fraction. For the 
one dimensional system of metallic sheets. 

FIG. 4. The photonic band structure for a square lattice of metal rods with a square cross-section. The filling fraction 
/ = 0.01, and w p d/2irc — 1. Calculated using the Order N method. 

FIG. 5. The photonic band structure for a square lattice of metal rods with a square cross-section. The filling fraction 
/ = 0.01, and w p d/2irc — 1. Calculated using the transfer matrix method. 

FIG. 6. The photonic band structure for a square lattice of metal rods with a square cross-section. The filling fraction 
/ = 0.04, and w p d/2irc — 1. Calculated using the Order N method. 

FIG. 7. The photonic band structure for a square lattice of metal rods with a square cross-section. The filling fraction 
/ = 0.04, and w p d/2irc — 1. Calculated using the transfer matrix method. 
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Figure 3 
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Figure 4 
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Figure 6 
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Figure 7 



